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Abstract 

Compliant environments can mediate interactions between me- 
chanically active cells like fibroblasts. Starting with a phenomeno- 
logical model for the behaviour of single cells, we use extensive Monte 
Carlo simulations to predict non-trivial structure formation for cell 
communities on soft elastic substrates as a function of elastic moduli, 
cell density, noise and cell position geometry. In general, we find a 
disordered structure as well as ordered string-like and ring-like struc- 
tures. The transition between ordered and disordered structures is 
controlled both by cell density and noise level, while the transition 
between string- and ring-like ordered structures is controlled by the 
Poisson ratio. Similar effects are observed in three dimensions. Our re- 
sults suggest that in regard to elastic effects, healthy connective tissue 
usually is in a macroscopically disordered state, but can be switched to 
a macroscopically ordered state by appropriate parameter variations, 
in a way that is reminiscent of wound contraction or diseased states 
like contracture. 



1 Introduction 



In order to develop and maintain tissues, cells in multicellular organisms 
have to interact with each other and the extracellular matrix (ECM). Cel- 
lular communication proceeds mainly through specific interactions provided 
by receptor-ligand binding. In solution, gradients in ligand concentration 
encode spatial information. For example, morphogen gradients guide cell 
differentiation P, and gradients of chemoattractants or chemorepellants di- 
rect cell motility (chemotaxis) 0. For cells adhering to each other or to the 
ECM, physical properties of the environment like topography and mechan- 
ics provide additional information supplementing specific biochemical cues. 
In particular, cells in their natural environment typically orient along fiber 
bundles of the ECM, a principle termed contact guidance [3 . In general, 
cells preferentially orient along directions with minimal curvature [3] . While 
contact guidance provides only a bidirectional cue for cell migration, unidi- 
rectional cues result from spatial gradients in adhesiveness (haptotaxis) [S], 
substrate rigidity (durotaxis) |Hj and substrate tension (tensotaxis) [7JIB]. 

Tissue cells like fibroblasts in the connective tissue constantly remodel 
their structural environment by degrading old and secreting new ECM. More- 
over they can exert forces large enough to actively reorganize the ECM after 
it has been laid down . Hence, cells not only respond to the physical prop- 
erties of their environment, but also actively modulate them. This results 
in indirect, matrix- mediated interactions between cells. In particular, cell 
traction-induced reorganization of collagen fibers can mediate a mechanical 
interaction between cells via contact guidance [PUl ITTj . In the same vein, 
elastic interactions between cells can result from stress and strain in the 
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ECM induced by cell traction. During recent years, the sophisticated use 
of elastic substrates has proven that cells indeed respond to purely elastic 
features in their environment, including rigidity, rigidity gradients and pre- 
strain [T21 Ed EH UH EES] • It now appears that many cell types, including 
fibroblasts, smooth muscle cells, and endothelial cells (but not neutrophils 
or neurons), respond to the mechanical properties of their environment with 
a common preference for large effective stiffness fTJEIJE]]. Here the term 
effective stiffness comprises both rigidity of and tensile strain in the environ- 
ment, which can be actively sensed by cells via mechanotransductory pro- 
cesses at cell-matrix adhesions [2111 12H 1221 I2H] • These mechanical cues play 
an important role in a variety of physiological processes, including develop- 
ment, tissue maintenance, angiogenesis, myotube fusion, wound healing and 
metastasis [2H 123 1213 12Z| • In particular, if coupled to cell division, stress and 
strain become major determinants of tissue morphogenesis [2E1I2I]]- In gen- 
eral, tissue function arises from the close relationship between cell behaviour 
in and material properties of the ECM [3U1 EH E2] ■ 

We have argued before that elastic interactions contribute to the way sin- 
gle cells position and orient themselves in compliant environments p£Hl34[ l35]. 
However, it has not been discussed in detail before how this individual be- 
haviour translates into collective behaviour of cell ensembles in compliant 
environments. In a recent short report we have shown that elastically inter- 
acting cells assemble into strings of cells which then leads to screening of the 
cellular traction patterns jSH]. Although this effect in principle suppresses 
macroscopic order, we have also found that macroscopic order can exist at 
high cell density, with an interesting competition between string- and ring- 
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like structures as a function of the Poisson ratio of the surrounding material. 
In our recent report, this effect has been especially analyzed for ordered (lat- 
tice) arrangements of cells, which can be decomposed into strings. Since 
cellular systems are intrinsically noisy, in this paper we focus on stochastic 
effects using extensive Monte Carlo simulations, which allow to predict col- 
lective effects in the presence of noise. In the following, we mainly consider 
the case of cells adhering to the top of a soft elastic substrate, because this 
setup might be the easiest way to experimentally verify our theoretical pre- 
dictions. Commonly used materials for elastic substrates are polyacrylamide 
or polydimethylsiloxane, which can be described by linear isotropic elasticity 
theory. Since we are interested in the interactions mediated by an elastic 
environment, we only consider situations without cell-cell contact, which ex- 
perimentally could trigger different cellular responses to mechanical signals. 
Therefore in our theoretical work we first fix cell positions and then relax cel- 
lular orientations using standard Monte Carlo techniques. Experimentally, 
this might be done by using microcontact printing, thus restricting cells to 
adhesive islands, or by using non-motile cell lines, which adhere at random 
positions and then do not move. By freezing in cell positions, we are also 
able to control cell density. For the position geometry, we consider ordered 
(lattice) arrangements, random perturbations around ordered positions and 
arrangements which are completely random. We calculate structural phase 
diagrams as a function of of elastic moduli, cell density, noise and cell po- 
sition geometry using Monte Carlo simulations. We then briefly investigate 
collective elastic effects in three dimensional elastic environments and finally 
discuss our results in the context of wound healing. 
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2 Model and simulations 

2.1 General concepts 

In order to sense the mechanical properties of their environment, tissue cells 
pull on it with actomyosin contractility. In many cases, this contractile 
mechanical activity is directed along the long axis of the cell body, espe- 
cially when confronted with a mechanically anisotropic environment. The 
mechanical action of such an anisotropic force pattern can be modeled as 
an anisotropic force contraction dipole Pij = Plilj, where the unit vector 
/ specifies cell orientation. In contrast to electric dipoles in electrostatics, 
which are vectors, force dipoles in elasticity theory are tensors of rank two. 
P specifies the dipole strength, which typically has \P\ = Fd « 10~ n J, 
corresponding to two opposing forces F = 200 nN separated by a distance 
d = 60 fim 137]. We will consider cells with anisotropic force patterns only 
and assume that the magnitude P is constant for all cells. A mechanically 
active cell generates an elastic strain field Uij(r), which again is a tensor 
of rank two. The strain induced by P can be calculated from the elastic 
Green tensor Gij(f,f') which describes the deformation of the medium at f 
caused by a point force at r*. In general Gij depends on material properties 
and boundary conditions. For translationally invariant situations, one has 
Gijif^r 1 ) = Gij(f — f 1 ). Elastic interactions between two cells result if the 
mechanical activity P^ of one cell responds to the elastic field Uij induced 
by another cell. In the framework of a general field theory elastic interaction 
can then be represented by a coupling of mj and P^. Since cells are active 
and moreover often show a regulated response, the exact form of the coupling 
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between and u y - cannot be predicted from first principles. In a first order 
approximation we may assume that P y - and linearly coupled. For me- 

chanically active cells like fibroblasts, experimental observations suggest that 
indeed they adopt positions and orientations in such a way as to effectively 
minimize the scalar quantity W = PijUij |34| 135] . For tensile strain, this im- 
plies that contractile cells actively align with the external field. Because they 
pull against the external stretch, cells reduce displacement. For fibroblast- 
like cells, this behaviour might have evolved in the context of wound healing, 
when cell traction is required to close wounds. A biophysical interpretation 
of the origin of the extremum principle regarding W is to note that W has 
dimensions of energy. Physically W can be interpreted as an extra energy 
contribution to the deformation energy required to build up a force dipole in 
the presence of strain tt y - |31] 135] . Using the analogy of a harmonic spring 
the basic idea can be explained easily: for a given spring constant K, it takes 
the work W = F 2 / (2K) to build up the force F. Therefore larger stiffness 
K corresponds to smaller work W required to reach the force F. In this 
way W may be interpreted as the inverse of an effective stiffness of the envi- 
ronment and the extremum principle in W corresponds to the experimental 
observation that cell-matrix contacts grow stronger in a stiff environment, 
thus eventually determining cell orientation. Using this extremum princi- 
ple, we have been able before to unify many diverse observations which have 
been reported for the organization of single cells on elastic substrates and 
in physiological hydrogels [3U|35]. F° r example, it predicts that single cells 
prefer to align in parallel and perpendicular to free and clamped surfaces of 
finite sized samples, respectively. It is however important to note that the 
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cellular response to the actively sensed effective stiffness of the environment 
is very different from the cellular response to cyclic external stretch with a 
1 Hz frequency, which is relevant for the physiology of lung tissue and the 
cardiovascular system. In this case, many cell types tend to orient away from 
the direction of stretch [3B1EH|, possibly to avoid the recurrent deformation 
of their cytoskeleton. 

For two cells, the above reasoning leads to the following effective interac- 
tion potential [3U EH] 

w = p l3 u l3 =Y1 p ^w-^i Gik(f ~ nPL (1) 

i,j,k,l 

where we consider only situations with translational invariance. The optimal 
cell configuration is described by the minimum of W in regard to cell posi- 
tioning and orientation. Since in general Gik scales as ~ l/(Er), where r is 
distance and E an elastic modulus, W scales as ~ P 2 / (Er 3 ) and has units of 
energy. Since in a linear material the cellular strain fields superimpose, the 
functional that describes elastic interactions of a system of N cells reads: 

1 N N 

"'^EE'"* ( 2 ) 

7=1 5^7 

where W^s is the interaction between two cells 7 and S as described in Eq. ^ 
The factor 1/2 is required to avoid double counting and the factor 1/N for per 
cell normalization. The optimal structure is given by the configuration which 
minimizes Eq. El as a function of all cellular orientations and positions and 
we may refer to this structure as the ground state, in analogy to interacting 
passive particles in physical systems. 
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2.2 Monte Carlo simulations 

In the presence of noise any system will deviate from its optimal state. In the 
cellular systems discussed here, noise results both from intracellular processes 
(like the intrinsic stochasticity of gene expression and signal transduction) 
and from heterogeneities in the material properties of the environment. In 
order to introduce a stochastic element into the structure formation process, 
we perform Monte Carlo simulations using the standard Metropolis algo- 
rithm to generate typical configurations in the presence of noise jlU]- This 
implies the use of an Gibbs ensemble, which in our case is the simplest choice 
given that we do not know the exact details of the cellular decision making 
process (in information theory, the Gibbs ensemble arises from maximizing 
the measure for disorder (Shannon entropy) under the constraint of a pre- 
scribed average energy). Starting from an arbitrary configuration, a cell is 
selected at random and its orientation randomly varied. The new configura- 
tion is always accepted when it decreases W. Otherwise it is accepted with 
the probability exp (— AW/ksT), where ks is the Boltzmann constant and 
T is temperature. In our context k B T represents a measure for the degree 
of stochasticity involved in cellular decision making. Other studies on mod- 
elling cellular structure formation have used the same concept before and 
found an effective value of k B T = 5 1CT 15 J [JTJ H21 EH| • Since this quantity 
comprises all molecular interactions between a cell and its environment, it is 
six orders of magnitudes larger than thermal energy. In our case, the com- 
petition between elastic and stochastic effects can be characterized by the 
reduced temperature 

^ k B TiiEtf 

T* = B p2 (3) 
7 



which measures the relative importance of noise with respect to the aver- 
age elastic interaction strength. Here b is the average distance between two 
cells which is related to the averaged cell density by (p) = 1/b 2 in two di- 
mensions (2D) and by (p) = 1/b 3 in three dimensions (3D). As we will see 
below, typical values in our simulations range from T* = 0.1 to 3.0. Using 
P = 10~ n J, E = 10 kPa, b = 100 pm and a typical value T* = 1, we obtain 
ksT = 3 10~ 15 J. Therefore our choice for the noise level is exactly in the 
same range as the earlier estimates ^T] H21 E3] • As the effective temperature 
is decreased towards cero, T* — > 0, W decreases and the ground states be- 
come increasingly favorable. This process is known as simulated annealing 
and is often used to identify optimal structures [44J . In the vicinity of an 
optimal structure, almost all Monte Carlo moves are rejected. In contrast, 
when T* — > oo, every Monte Carlo move is accepted and the structures get 
disordered. One Monte Carlo sweep corresponds to N such Monte Carlo 
moves. After thermal equilibration (typically about 10 3 — 10 4 Monte Carlo 
sweeps), the Metropolis algorithm samples the important configurations typ- 
ical for a given temperature T*. One can then calculate statistical properties 
of structures at defined noise level by averaging the quantity of interest over 
many configurations generated by the Metropolis algorithm. 

To study stochastic effects using our Monte Carlo simulations we typically 
consider iV « 1000 dipoles. In order to minimize the effects of boundaries 
and finite size, we apply periodic boundary conditions (pbcs) , such that each 
dipole has the same number of next neighbors and experiences the same 
local geometry. We implement pbcs using the minimal image convention, 
i.e. we only consider the interactions of the dipole with its N — 1 nearest 



8 



(image) dipoles [30]. Note that in the two-dimensional situation of cells 
on top of an elastic half space, the 1/r 3 interaction effectively constitutes a 
short-ranged potential, because its area integral does not diverge with system 
size (f L rdr(l/r 3 ) ~ 1/L where L is system size). Thus, boundary effects 
are expected to play only a minor role in this case and the minimal image 
convention is a good approximation. 

2.3 Additional assumptions 

We now make some additional assumptions which will simplify our subse- 
quent calculations and which will allow for a direct comparison of our the- 
oretical calculations with appropriate experiments. First, we assume that 
the environment can be described by isotropic linear elasticity theory, which 
is a reasonable assumption for the synthetic substrates like the ones made 
from polyacrylamide or polydimethylsiloxane commonly used to study me- 
chanical effects in cell culture |45| I37j. Thus there are two elastic constants: 
the Young modulus E describes the rigidity of the material and the Poisson 
ratio v the relative importance of compression and shear. In particular we 
mainly consider the situation of cells adhering to the top surface of a thick 
elastic film. Such a situation can be theoretically represented by an elastic 
half-space geometry with dipolar orientations constrained to the x-y-plane. 
Therefore the Young modulus E and the Poisson ratio v are 3D quantities. 
The maximal value for v is 1/2, e.g. for strongly hydrated polymer gels. If 
such a material is tensed in one direction, the shear mode is excited and it 
contracts in the perpendicular directions (Poisson effect). For common ma- 
terials, the minimal value for the Poisson ratio is v = 0, e.g. in dehydrated 
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fibrous polymer gels. Then the volume mode prevails and uniaxial tension 
does not translate into lateral contraction. For cells exerting tangential forces 
on top of an elastic substrate one can use the Boussinesq Green function of 
an elastic halfspace 46: to specify Eq. [T]to [35J: 



where r is the distance between cells and the cellular orientations 9, 6' and a 
are defined via the scalar products cos 9 = l-f, cos 9' = V • r and cos a = l-V . 
f is given by 



where a\ = z/(l + v)/ (nE) and a<i = (1 — v)jv. 

Secondly, because here we focus on elastic effects, we want to avoid cell- 
cell contacts, which are known to change the mechanical state of adhering 
cells [TSj. Therefore we attribute an exclusion disc of radius a to each cell. 
Moreover in most of our simulations cell positions are frozen and only ori- 
entational degrees of freedom are considered. Experimentally, this situation 
might be achieved by using non-motile cell lines or appropriate drugs to sup- 
press cell motility. Alternatively, one might use microcontact printing to 
prepare well-defined adhesive islands constraining cell positions. 

Finally, while the focus of our paper is on structure formation on planar 
synthetic substrates to allow a direct comparison of theory and experiment, 
we also briefly consider cells in a three-dimensional (3D) environment. For 
computational simplicity we keep the assumption of linear isotropic. For cells 




(4) 
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embedded in an infinite 3D isotropic elastic environment, the interaction law 
stays the same as Eq.0 but with different constants a% = (1+u) / (SttE(1— u)) 
and <i2 = (3 — 4z/) which now follow from the Kelvin solution for the full 
elastic space |46j . In 3D, the elastic interactions are truly long-ranged, in the 
sense that now the volume integral over the elastic interaction diverges with 
system size (J L r 2 dr(l/r 3 ) ~ InL where L is system size), thus boundary 
effects become more important. 

3 Results and Discussion 

3.1 Structure formation for dipolar lattices 
3.1.1 Phase diagrams under thermal noise 

We first consider structure formation in the presence of noise when cells 
adopt regular positions on an infinite square (s) and hexagonal (h) lattice. 
In general we find a strong dependence of structure formation on lattice ge- 
ometry, Poisson ratio v and noise level T*. In particular our simulations 
show an interesting competition between string-like and ring-like structures 
at low noise levels, which we have also found before in a detailed analysis 
of optimal structures [UHl- For the square lattice low noise structures are 
always string-like. The lattice geometry favors the formation of strings di- 
rected along one of the principal lattice vectors n = (1, 0) or n = (0, 1) for 
all v. The limit v — > shows structural degeneracy, that is, additional struc- 
tures, e.g. parallel strings aligned along n = (1, 1), become equally favored as 
parallel (l,0)-string structure. Enhanced structural degeneracy of string-like 
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structures in the limit of vanishing Poisson ratio occurs also on the hexag- 
onal lattice. Most interestingly, however, on the hexagonal lattice string- 
and ring-like structures exchange stability as a function of Poisson ratio. For 
v > 0.32 the Poisson effect favors ring-like structures over strings, with rings 
composed of four dipoles each rotated by 90 degrees relatively to its neighbor. 
On both lattice types the orientational order imposed by elastic interactions 
is gradually destroyed with increasing noise levels. When elastic interactions 
dominate, dipoles typically fluctuate weakly around their optimal orienta- 
tions forming long-ranged ordered structures. At intermediate noise levels, 
long-ranged order is decreased, but there are still ordered domains which 
locally adjust to elastic signals. As noise increased further, domain size de- 
creases and local fluctuations become stronger, eventually giving eventually 
rise to completely random dipole arrangements. 

The loss of long-ranged order with increased noise can be quantified by 
calculating the temperature behavior of appropriate order parameters. For 
string-like structures all dipoles point along a common direction n. A suitable 
parameter to quantify the degree of global alignment therefore seems to be 
the angle f3 between dipole orientation / and n, where cos (3 — l-n. Note that 
force dipoles have a bipolar symmetry and thus dipole orientations / and —I 
are equivalent. Thus, the average (cos (3) = 0. However, the average (cos 2 /3) 
becomes one when all dipoles point along n, while for an isotropic structure 
(cos 2 f3) = \. Thus, we define an alignment order parameter p as 

p = 2(cos 2 /3) - 1, (6) 

which is non-zero only for globally aligned structures and zero otherwise. 
From our Monte Carlo simulations we compute p as a function of T* and 
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Poisson ratio for both lattice types. As shown in Fig.^i, on the square lattice 
a globally aligned structure p — > 1 always develops regardless of Poisson ratio 
when the noise level becomes sufficiently low. The transition temperature T* 
from a state without into a state with long-ranged orientation correlations, 
i.e. where p 7^ for the first time, decreases approximately linearly with v. 
This might be attributed to the general enhancement of elastic interactions 
with increasing Poisson ratio W ~ (1 + v). In Fig. we show how the 
probability distribution to find a dipole with a given orientation f3 evolves 
with increasing noise levels. At low T* almost all dipoles point along n 
and the orientation distribution is strongly peaked; with increasing T* the 
distribution broadens and finally for T* > T* all orientations are - on a 
global scale - equally likely (circle). In contrast to the square lattice, where 
long-ranged order always implies alignment, the type of long-ranged order 
developed on a hexagonal lattice depends on Poisson ratio. Fig. shows 
that an aligned structure develops only for v < u c « 0.32. Interestingly, in 
this case the relation between the transition temperature T* and Poisson ratio 
is also reversed compared to the square lattice: T* stays roughly constant 
for v < 0.2 and then drops quickly around v « v c . Moreover, towards v c 
the transition seems to become more discontinuous and in our simulations 
p suddenly jumps from a low to a large value, indicating bistability and a 
first order phase transition. The T* scaling with v might be explained by 
simultaneous overall decrease in elastic signal strength with decreasing v 
within strings compensated by an increase in lateral string-string coupling 
favoring global alignment at small v. For larger Poisson ratio {y > 0.4) we 
find p = at all temperatures. At low T* < T* a macroscopic isotropic 
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long-ranged ring-like order exists build up by small rings composed of four 
dipoles oriented at 90 degrees with respect to each other. That is, at T* < T c 
there are two equally large subpopulations of dipoles oriented at 90 degrees 
with respect to each other. This corresponds to the two peaks in the radial 
orientation distribution shown in Fig.^i (y — 0.5). For v = 0.5 the transition 
temperature into the ordered ring-like structure T* is slightly larger than 1 
and significantly lower than T* rs 2 observed for the string structures on the 
square lattice at equal dipole densities. 

In El and |3] we summarize our results for ordered structures in the form 
of v — T* phase diagrams for square and hexagonal lattices, respectively. 
We have included representative snapshots of typical dipole configurations 
at noise levels well below T* in the globally ordered state, approximately 
at the phase transition T* and well above T* in the disordered phase for 
v = 0.1 and v = 0.5, respectively. It is important to note that despite the 
loss of global order at elevated noise levels T* > T*, elastic interactions still 
locally manifest themselves in characteristic domain formation. For example, 
at T* = 3 elastic interactions still cause local order on a square lattice and 
one finds short string-like domains preferentially oriented along the principal 
lattice vectors. On incompressible substrates {y = 0.5) we furthermore often 
observe cooperative excitations and the formation of ring-like domains. This 
may reflect the competition between ring- and string-like domains which is 
modulated by the Poisson ratio. In fact, as we have shown before, at v = 0.5 
ring-like structures are only slightly disfavoured with respect to the string-like 
structures jHH]- Nevertheless, it may seem surprising that ring-like domains 
are excited so easily already at very low T* even before coexistence of the 
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equally favorable (1,0) and (0,1) domains occurs (see snapshot at T* = 1). 
The reason is probably a large domain interface penalty for two (1,0) and 
(0,1) string domains compared to the interface penalty of a string and an 
adjacent ring-domain. 

For the hexagonal lattice local structural characteristics of elastic interac- 
tions under noise are again string-like domains preferentially oriented along 
directions specified by lattice geometry. For large u, a characteristic signa- 
ture for elastic interactions is the formation of 4-cell-ring or ladder-domains 
with cells having approximately perpendicular orientations with respect to 
each other. Thus, although experimentally noise might be too strong to allow 
for global ordering, elastic interactions might be still detected by their local 
signatures of characteristic domain formation. 

3.1.2 Effect of weak positional disorder 

Experimentally cells may not adopt perfect lattice positions, in particular 
when adhesive islands are used which are large compared to cell size in order 
to minimize the effects of island shape on cellular force distribution jl7| HH] • 
We therefore investigate the effect of positional fluctuations at low thermal 
noise intensities. For this purpose we randomly displace the dipole positions 
within a circle of radius r around the lattice positions before relaxing the 
orientations. The degree of positional disorder can be quantified by the ratio 
/ = r/b of the radius with respect to the lattice constant b. In the following 
we focus on the square lattice and v = 0.5, since synthetic polymer gels 
typically have v m 0.5. For this value of the Poisson ratio, Fig. ^ predicts a 
lattice geometry dependent transition from a string-like structure on a square 
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lattice to a ring-like structure on the hexagonal lattice. However, rings are 
thermally excited easily on the square lattice. Hence the string-like structure 
might also be easily destabilized under positional disorder. 

In Fig. 0] we show Monte Carlo snapshots at a low reduced temperature 
(T* = 0.1) for / = 0.1,0.25,0.5. The simulations show that a 10% deviation 
from the lattice positions has only a minor effect on orientational ordering. 
For / = 0.25, string-like as well as ring-like domains form and long-ranged 
orientational order disappears. As for thermal noise, ring-domains typically 
localize to interfaces between string-like domains with perpendicular orienta- 
tions. With increasing / domains shrink and structures appear increasingly 
disordered. Thus, moderate deviations from the square lattice positions are 
not sufficient to destabilize string-like structures with respect to ring-like 
structures on incompressible substrates. Possibly because under uniform po- 
sitional noise dipoles still find themselves in a local fourfold coordination 
with other dipoles (that is, each dipole has typically for next neighbors as 
on the square lattice), string domains are stabilized. We conclude that posi- 
tional disorder affects structure formation in a similar way as increasing the 
temperature T* on a perfect lattice. 

3.2 Phase behavior on homogeneous substrates 

On a homogeneous substrate, non-motile cells usually adhere more or less at 
random positions. Therefore we next study typical structures on elastic sub- 
strates with completely randomized but fixed positions as described above. 
Cellular structure formation on elastic substrates is now governed mainly by 
three control parameters, namely reduced temperature T*, Poisson ratio v 
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and cell density. Since in our model each cell is characterized by an exclusion 
radius a, we introduce a reduced density 

P = -JjT, (7) 

which is a dimensionless variable describing the ratio of the area occupied by 
N circular disks of radius a to the area of the (simulation) box with length 
L. 

Fig. El shows typical snapshots of structures at T* = 0.1 for cells on an 
elastic substrate with v = 0, 0.25, 0.35, 0.5 (top-bottom) for p* = 0, 0.4, 0.5 
(left -right). At low densities cells predominantly optimize locally the in- 
teraction between them by forming short strings with no long-range corre- 
lation between string-like clusters. We attribute this finding to the strong 
tendency of dipoles to form strings and the strong elastic screening of string- 
string interactions j^S]. This results in rather robust pattern formation that 
does not differ qualitatively as the Poisson ratio is varied, see Fig. ISk.. One 
expects that these patterns represent typical cellular structures formed by 
strongly interacting cells when cells in dilute concentrations are suspended 
on an elastic substrate and adhere at random positions (p* — > 0). With 
increasing p* the respective structures at low noise intensity show a strong 
dependence on the Poisson ratio v. For incompressible substrates [y = 0.5) 
we find isotropic ring-like structures, see Fig. EpVc). With decreasing Pois- 
son ratio string-like patterns emerge. For v = 0.35 we find coexistence of 
string-like and ring-like domains, compare Fig. Exilic). For substrates with 
Poisson ratio below v < 0.32 string-like structures dominate at intermediate 
densities. With increasing density strings start to interact and domains of 
aligned parallel strings form which increase in size with increasing p*, see 
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Fig. |^IIb,IIc). Highly compressible substrates [y — > 0) favor cell alignment 
along a common direction at high cell densities and in Fig. EJdc) we find a 
globally aligned structure. Here, the rotational symmetry of the structure 
is spontaneously broken along an arbitrary direction in space despite lack- 
ing long-ranged positional order. In this respect we may speak of a nematic 
structure in analogy to nematic phases formed by liquid crystals which are 
characterized by orientational order but positional disorder. 

As before one can quantify the development of anisosotropic long-ranged 
order with the order parameter p defined in Eq. In contrast to the lat- 
tice structures where preferred string direction n was determined by lattice 
geometry, the orientation of the director n on homogeneous substrates is ar- 
bitrary. This is taken into account by defining a two-dimensional analog of 
the nematic order parameter p used in the theory of liquid crystals jlH] • We 
first introduce the ordering matrix Qij = ^$^a=i(^i — §<%)> where l a is 
the orientation vector of the a'th particle and the sum runs over all particles 
in the simulation box. The largest eigenvalue A of the symmetric ordering 
matrix Q corresponds then to the order parameter p = 2A. As before p 
measures the degree of orientational order with respect to the current di- 
rector n, which is the corresponding eigenvector to the maximal eigenvalue. 
p = 1 implies nematic and p = isotropic structures. Since we freeze in 
positions, we first thermally average p for a fixed configuration of the dipole 
positions and subsequently average over at least 20 random position config- 
urations obtained for the same p*. In Fig. |H]we show the quantitative effects 
of the reduced density p*, Poisson ratio v and reduced temperature T* on 
nematic ordering. At constant temperature T* = 0.1 a nematic structure 
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forms beyond a critical density p*(v) and below a critical value of the Pois- 
son ratio, see Fig. The degree of structural alignment increases with 
increasing p* and approaches p — *■ 1 toward the maximal p* = 0.907, which 
is the maximal packing density achieved for disks arranged on a hexagonal 
lattice. For v > 0.4 no nematic ordering p = occurs at any density p* or 
temperature T*, because ordered structures are ring-like. In addition to the 
density-dependent ordering transition, the nematic structure is also destabi- 
lized by increased noise. In Fig. \§jp we plot p for v = at different values 
of the reduced temperature T*. The critical density p* to obtain a nematic 
structure increases with T* and p = for any any p* when T* exceeds a 
critical T* ~ 1, which corresponds approximately to T* measured on the 
hexagonal lattice. 

In Fig. we show the phase diagram in the z/-p*-plane for T* = 0.1. 
Diamonds mark p < 0.4, squares p > 0.4 and the dashed line denotes our es- 
timate for the isoline p = 0.4. We find three different phases: a high-density 
anisotropic string-like structure, a high-density isotropic ring-like structure 
and a low density disordered structure. The general features of this phase 
diagram are very similar to the thermal phase diagram for hexagonal lattices 
shown in Fig. El but with the role of T* and (inverse) density p* exchanged. 
Note, that at p* = 0.907 dipoles in fact adopt hexagonal lattice positions, 
because disks must order into a hexagonal lattice due to the non-overlap 
constraint. Thus, reducing the density from p* = 0.907 in our model is in 
some sense equivalent to introducing position fluctuations around hexagonal 
lattice positions (see above), which qualitatively has similar consequences as 
raising the effective temperature. These results also show the importance 
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of local position geometry and dipole coordination. At high p* dipoles have 
six next neighbors and all interactions are equal in magnitude. Because 
elastic interactions are strongly anisotropic and dependent on Poisson ratio, 
non-trivial structure formation as a function of local dipole coordination and 
Poisson ratio results. Decreasing p* the number of next neighbors decreases 
and at low p* dipoles typically only have one or two next neighbors, which 
favors dipole alignment and string formation. Due to screening of the trac- 
tion patterns in a string, adjacent strings hardly interact and we obtain a 
largely disordered structure even at low T*. In our calculations the degree of 
position correlation is controlled by p*, which therefore plays a similar role 
as lattice geometry before. In contrast to relative density effects described 
by p*, absolute density effects can be subsumed into the effective tempera- 
ture T*, see Eq. El This is due to the fact that by increasing cell density 
the screening efficiency of string-string interactions also increases. The only 
way to overcome screening and introduce lateral coupling between dipoles is 
to insure that dipole-dipole distances within a string and the distance to an 
adjacent string are comparable. Increased dipole coordination, however, is 
controlled by p* rather than by p. 

We note that the order- disorder transition with p* might also represent 
a spin-glass transition. Such transitions are common for magnets, where 
spins represent magnetic moments. In magnets global ordering of spins can 
result in a ferromagnetic state, where all spins point along a common direc- 
tion and a macroscopic magnetic field builds up. Anti-ferromagnetic order 
results when spins are globally ordered but with equal numbers of up and 
down spins, such that no macroscopic polarization builds up. A spin-glass 
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transition occurs when e.g. due to disordered positions of spins no ferro- 
or anti-ferromagnetic order develops even at low temperatures. In analogy, 
dipoles in our model represent spins. The aligned structure therefore might 
be called a ferroelastic phase, because a macroscopically anisotropic polariza- 
tion stress develops spontaneously. In analogy, the ring-like structure might 
be called an antiferroelastic phase, because the macroscopic stress exerted by 
the dipoles remains isotropic. When p* < p c , no global order develops even 
in the limit T* — > 0; instead, many equally favorable states exist. Because 
we do not allow dipoles to adjust positions, they cannot relax into an ordered 
ferro- or anti-ferroelastic state. Thus, the disordered low temperature phase 
might also represent a spin glass. There is an interesting difference between 
thermally disordered states and density disordered states at low T. In both 
cases, dipoles can be aligned by application of an external field. The former 
might be called a paraelastic phase because it relaxes into the unpolarized 
disordered state after the external field is switched off with a single time-scale 
t oc (T — T c ) 7 . In contrast, the later might represent a spin-glass state and 
relaxation into the disordered state cannot be described by a single power 
law. In fact a characteristic feature of spin glasses is a fast relaxation to a 
remanent polarization and then a slow relaxation into the unpolarized state, 
which may exceed experimental timescales. 

4 Elastic effect in 3D environments 

Although cell behaviour in 3D is not the focus of this paper, we now spell 
out a few predictions resulting from our model in this case. Calculations 
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for cells positioned on a simple cubic lattice suggest that in incompressible 
environments {y = 1/2) the optimal structure is effectively isotropic with a 
hedgehog-like unit cell, where all dipoles at the corners point to the cube's 
center, see Fig. while for v — spontaneous symmetry breaking along a 
principal lattice lattice vector occurs, see Fig. [7)3. This indicates that sim- 
ilar transitions between ring-like and string-like structures as a function of 
Poisson ratio v exist in 3D as they do in 2D. String-like structures might 
also be favored by anisotropies in the effective mechanical properties of the 
environment. For example, fiber alignment or external strain fields might 
cause cell alignment with the external perturbation, which could be further 
stabilized by elastic interactions between cells. In Fig. UJ: we show a corre- 
sponding snapshot of a Monte Carlo simulation with 100 cells modeled as 
hard spheres with an elastic dipole moment at their center. Here we allowed 
for both orientational and positional degrees of freedom (T* = 2, v — 1/2), 
where the position dependence of W from Eq. |2]is used for additional Monte 
Carlo moves. Moreover we applied a homogeneous strain field along the 
^-direction. Cells respond by an alignment with the external field and the 
formation of strings as has been observed experimentally [SU] • Without exter- 
nal field but with mobile cells we typically observe formation of "networks" 
of cells at low temperature and intermediate densities and disconnected mi- 
correlated strings at elevated noise levels and low densities, respectively. In 
general, we expect that the low to intermediate density regime of a three 3D 
phase diagram for mobile elastically interacting force dipoles will resemble 
the phase behavior of electric dipolar fluids I54j . In both systems string 
formation dominates and string-string interactions are screened in the same 
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way [55]. At high dipole densities our lattice calculations in 3D suggest that 
phase transitions as a function of Poisson ratio may occur and we find evi- 
dence that again a low Poisson ratio of the material tends to favor global cell 
alignment. 

In finite 3D environments the geometry and boundary condition of the 
sample may also affect structure formation. The presence of a boundary mod- 
ifies the direct elastic interaction between cells by boundary induced strain 
fields, which, depending on boundary condition, introduce either attractive 
or repulsive contributions to the elastic interaction plus a direct interaction 
term with the boundary |33] ■ An instructive example for the competition 
between the direct interaction of a cell with the boundary and cellular in- 
teractions is the elastic half space with a clamped boundary. When cells are 
located close to the boundary, the direct interaction with a clamped bound- 
ary favors cellular orientations pointing toward the surface [HUES]- O n the 
other hand, interactions between cells favor the formation of strings and thus 
parallel orientations. Our calculations suggest that the transition between 
these two configurations is a function of the ratio a = b/d, where d is the 
distance to the boundary and b is the distance between cells, and the number 
of interacting cells N. Above a critical value a c ~ 2.3 the direct interaction 
with the boundary dominates and cells are expected to point toward the 
surface. For a < a c orientations in parallel to the clamped surface are fa- 
vored when sufficiently many cells are present, N > N c (a). In Fig. [7|i we 
plot N c as a function of a. Thus, collective effects may alter preferred cell 
organization close to clamped boundaries. In contrast for a free surface the 
direct interaction with boundary favors parallel alignment, which is further 
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stabilized by collective elastic effects between cells. This may explain the ro- 
bust parallel alignment of cells with respect to free surfaces which has been 
observed experimentally |oT| I52|. 

5 Conclusion and Outlook 

Using a simple mathematical model and extensive Monte Carlo simulations, 
we have analyzed structure formation due to elastic effects for fibroblast-like 
cells with anisotropic force patterns. Our model predicts in a quantitative 
way how structure formation is controlled on soft elastic substrates by the 
Poisson ratio u, the reduced cell density p*, the relative strength of elastic 
interactions with respect to noise (represented by the reduced temperature 
T*) and the geometry of cell positioning. Up to intermediate cell densities, 
we predict the formation of uncorrelated strings of cells (paraelastic phase). 
This finding can be explained by noting that the cellular traction patterns 
in strings screen each other, as we have demonstrated before by an analyt- 
ical calculation [SB]. For high cell densities, we find an interesting compe- 
tition between isotropic ring-like and anisotropic string-like structures. The 
string-like (ferroelastic) phase is reminiscent of the nematic phase for liquid 
crystals. In contrast, the nematic order parameter vanishes in the ring-like 
(anti-ferroelastic) phase, because local ordering principles preclude nematic 
order on the macroscopic scale. 

The Poisson ratio v allows to switch between these phases. It strongly 
affects cellular structure formation due to the non-trivial way in which stress 
and strain propagates in the elastic environment. The antiferroelastic phase 
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is expected to occur on synthetic elastic substrates made from e.g. polyacry- 
lamide or polydimethylsiloxane, which are characterized by v ~ 1/2. In order 
to obtain the nematic phase, new kinds of polymer gels for cell culture are 
needed, with considerable lower values for the Poisson ratio (namely below 
v c = 0.32). One possible route might be the design of biocompatible polymer 
gels with large meshsizes but little hydration, thus avoiding the incompress- 
ibility effect of bound water. As it is an open issue on which time and length 
scales cells sense the mechanical properties of their environment [19 j , design 
and testing of new materials is crucial for progress in this field. 

Like the Poisson ratio, also the other determinants of cellular structure 
formation are in principle accessible experimentally. The Young modulus E 
directly affects the reduced temperature and can easily be varied in experi- 
ments. Nevertheless it does not appear to be a reasonable control variable 
for experiments, mainly because the predicted structure formation might be 
expected to take place only in a small range of .E-values, namely the ones 
corresponding to physiological rigidities (around kPa). Other control param- 
eters of interest are cellular contractility (which affects T* through the force 
dipole moment P and which might be varied for example by administering 
LPA, which stimulates Rho-mediated contractility) and the average distance 
b between cells (although, in order for the force dipole approximation to re- 
main valid, the typical distance between cells should not be smaller than 
the typical cell length). The most important experimental control param- 
eter next to the Poisson ratio v might be the geometry of cell positioning, 
which can be controlled with patterning techniques (including microcontact 
printing in two dimensions). 
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Elastic effects are likely to contribute to large scale tissue organization 
and to the structure-function relationship of tissues. Our results indicate 
that for sparsely populated tissues, elastic interactions, despite their long- 
ranged character, might have only a rather local effect. This is a result of 
the local tendency for cells to align into strings combined with the resulting 
screening of elastic fields. In this regime the system remains macroscopically 
disordered over a wide range of variations in material properties and cell 
density, effectively making the composite material of matrix and cells robust 
against perturbations. It is well known that the sparsely populated connec- 
tive tissue is a rather disordered structure. Our simulations suggest that 
without macroscopic external fields elastic interactions between cells favor 
isotropic disordered structures. In this case the average stress in the tissue 
resulting from cells is isotropic, that is aij = PpSij, where p is the cell density. 
The situation changes when an external field is applied (e.g. resulting from a 
wound) . Then highly aligned structures may prevail. In this case cells orient 
their mechanical activity in such a way that the average cellular stress in the 
medium = p(Pij) becomes anisotropic in response to homogeneous uni- 
axial stress, that is, the cellular stress is directed opposite to the externally 
applied stress (e.g. as to close the wound). It is an interesting question what 
happens after the external alignment field has decreased to zero. Our cal- 
culations suggest that elastic interactions may stabilize an aligned structure 
even in the absence of external fields. In particular, disordered tissues (like 
the connective tissue) may become aligned due to an external field. Due to 
elastic interactions the transition from the aligned state back into the normal 
disordered state upon switching off the external field could be dramatically 
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slowed down, in analogy to similar effects observed in magnetic spin glasses. 
This could have implications for wound healing and scar formation. Polariza- 
tion also occurs spontaneously in our model. Interestingly, this state reminds 
of the disease state called contracture, when certain tissues (like skin or mus- 
cle) start to tighten up, eventually preventing movement of body parts. It is 
an interesting speculation whether the mechanical properties of extracellular 
matrices are altered by (diseased) cells in such a way as to effectively induce 
variations in the Poisson ratio, which could trigger transitions between ferro- 
and anti-ferroelastic states. In order to investigate this point in more detail, 
our model should be extended to include viscoelastic and plastic effects as 
well as fiber degrees of freedom. 

Mechanotaxis also plays an important role for dynamic rearrangement of 
groups of cells in tissues, including development, wound healing, capillary 
sprouting and tumor growth. The work presented here builds on experimen- 
tal observations which suggest that cells migrate toward high strain areas 
and orient their mechanical activity in such a way as to pull back in response 
to external tensile strain. Such tensile strain is certainly present e.g. close 
to wounded areas or tumors. Future studies will show in which sense our 
model can be adapted to these more specific situations. The current study 
shows in a general way that many interesting structure formation processes 
arise from mechanical effects in communities of mechanosensitive cells and 
therefore might help not only to understand the corresponding physiological 
situations, but also to program new tissue functions by guiding the design of 
novel material with appropriate mechanical properties. 
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Figure Captions 



Figure 1: Effect of noise on global alignment order, (a) The order pa- 
rameter p for alignment on a square lattice as a function of noise T* for 
Poisson ratio v = 0,0.1,0.2,0.3,0.4,0.5 (from left to right). For low noise, 
aligned structures develop for all values of the Poisson ratio, (b) The radial 
distribution of the dipolar orientations on a square lattice with T* = 1.5, 2, 3 
for Poisson ratio v — 1/2. The two peaks correspond to an aligned structure 
with bipolar symmetry below T£ = 3. (c) On a hexagonal lattice, p = 
for v = 0.4,0.5 and p — > 1 for v = 0.3,0.2,0.1,0 from left to right. At low 
noise, aligned structures are unstable for large values of the Poisson ratio, 
(d) The radial distribution function on a hexagonal lattice with T* = 0.5, 1, 2 
for Poisson ratio v — 1/2. The four peaks correspond to a ring-like structure 
below T* = 2. 

Figure 2: Schematic phase diagram as a function of noise level T* and Poisson 
ratio v for the square lattice. For illustration, typical Monte Carlo snapshots 
are added. Both decreasing T* or increasing v favors the ordered string-like 
phase. For v — 0.1, the dipoles typically weakly fluctuate around the string- 
like ground state. For v = 0.5, the dominant fluctuations around the optimal 
state are ring-like. 

Figure 3: Schematic phase diagram as a function of noise level T* and Pois- 
son ratio v for the hexagonal lattice. For illustration, typical Monte Carlo 
snapshots are added. Ordered phases are favored by decreasing T*, but the 
ordered phases are string- and ring-like for small and large values of z/, re- 
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spectively. 



Figure 4: Snapshots of Monte Carlo simulations at reduced temperature 
T* = 0.1 for cellular force dipoles randomly displaced from an ideal square 
lattice for / = 0.1,0.25,0.5 from left to right for v = 0.5 (a-c) and v — 
(d-f), respectively. Increased positional disorder / destroys long-ranged ori- 
entational order in a similar way as an increase in T* on a perfect lattice. 

Figure 5: Snapshots of Monte Carlo simulations at reduced temperature 
T* — 0.1 for N = 1024 cellular force dipoles at random, but fixed positions 
on an elastic substrate with v = 0,0.25,0.35,0.5, respectively (top-bottom). 
The reduced area density increases from left to right as p* = 0, 0.4, 0.6, re- 
spectively, while the average cell density (p) remains constant. 

Figure 6: Averaged nematic order parameter p as a function of model param- 
eters, (a) Simulation results for v — 0, 0.25, 0.3 (left-right) and v = 0.4, 0.5 
(bottom), respectively. Thus nematic ordering occurs only on compressible 
substrates and above a critical density, (b) Simulation results for v = at 
T* = 0.1,0.6,1.1 (top-bottom). Thus nematic ordering disappears above a 
certain noise level, (c) Phase diagram for T* = 0.1 for dipoles on elastic 
substrates obtained by Monte Carlo simulations. All points below diamonds 
have p < 0.4 and all points above squares yield p > 0.4. The long dashed 
line is our estimate for the isoline with p = 0.4. The horizontal dashed line 
marks the maximal p* possible for a hexagonal lattice. 
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Figure 7: Collective effects in 3D elastic media. (a,b) The optimal struc- 
ture in a cubic lattice depends on the Poisson ratio v. In incompressible 
media [y = 1/2) the isotropic hedgehog structure is most favorable (a) while 
in highly compressible media [y = 0) elastic interactions favor aligned struc- 
tures (b). (c) Cells in external strain fields form strings running along the 
direction of stretch due to interactions with external strain and elastic cell- 
cell interactions, (d) Collective effects can modify preferred cell organization 
close to a clamped boundary, which for single cells is perpendicular. For 
small intercellular distances (a = b/d < 2.3) and sufficiently many cells in a 
string (N > N c ), the parallel orientation becomes more favorable. 
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Figure 1: Bischofs and Schwarz 
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Figure 2: Biscr^s and Schwarz 
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Figure 3: Bischofs and Schwarz 
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Figure 6: Bischofs and Schwarz 
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Figure 7: Bischofs and Schwarz 
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